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A first principles phase diagram calculation, that included van der Waals interactions, was performed for the 
3D bulk system (1 — X) • M 0 S 2 — {X) • MoTe 2 . Surprisingly, the predicted phase diagram has at least two 
ordered phases, at X ^ 0.46, even though all calculated formation energies are positive; in a ground-state 
analysis that examined all configurations with 16 or fewer anion sites. The lower-temperature /-phase is 
predicted to transform to a higher-temperature /'-phase at T ^ 500//, and /' disorders at T ^ 730//. Both 
these transitions are predicted to be first-order, and there are broad two-phase fields on both sides of the 
ordered regions. Both the /- and /'-phases are predicted to be incommensurate i.e. aperiodic: /-phase in 
three dimensions; and /'-phase in two dimensions. 


I. INTRODUCTION 


Recently there has been great interest in two- 
dimensional (2D) transition metal dichalcogenide (TMD) 
materials such as M 0 S 2 , MoSe 2 and MoTe 2 , their 
solid solutions, and related 2D materials^’^. Tradition¬ 
ally, M 0 S 2 has been used as a dry lubricant^ that 
is stable up to 623 K. Currently, interest is focused 
on applications as: band-gap engineering materials^’^; 
nano-electronic devices^’^”^; photovoltaic devices^’^^; val- 
leytronics applications^^’^^; 2D building blocks for elec¬ 
tronic heterostructures^^; and as sensor materials^^. 

The individual, three-atom-thick, 2D-layers of the bulk 
system are bonded by van der Waals forces, hence these 
forces influence bulk and multilayer synthesis and there¬ 
fore anion order-disorder and/or phase separation in solid 
solutions. The results presented below, for 3D bulk 
M 0 S 2 — MoTe 2 , imply that van der Waals interactions 
may strongly affect phase stabilities, either between ad¬ 
jacent layers in bulk or few-layer samples, or between 
monolayers and heterogeneous substrates. 

Of the three quasibinary solid solutions {M 0 S 2 — 
MoSe 2 ^ MoSe 2 — MoTe 2 , M 0 S 2 — MoTe 2 ) (1 — X) • 
M 0 S 2 ~ {X) ' MoTe 2 has the greatest difference in an¬ 
ionic radii (R5'=1.84 A; RTe=2.21 A)^^, which suggests 
that it is the most likely to exhibit interesting solution be¬ 
havior. One expects a simple miscibility gap as reported 
by Kang et al.^ for monolayer M 0 S 2 — MoTe 2 , hence 
the prediction of two configurational entropy {Scon) sta¬ 
bilized incommensurate, i.e. aperiodic, phases is ex¬ 
traordinary (stable phases that have positive formation 
energies must be entropy stabilized). 
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II. METHODOLOGY 
A. Total Energy Calculations 

Total structure energies, AEstr were calculated for 
fully relaxed M 0 S 2 , MoTe 2 (2H-structure, space group 
PQs/mmc, AB-stacking of three-atom-thick layers), and 
for 233 Mo^+n(SmTen )2 supercells. The Vienna 
ab initio simulation program (VASP, version 5.3.3 ^^47^ 
was used for all density-functional theory (DFT) calcula¬ 
tions, with projector augmented waves (PAW) and a gen¬ 
eralized gradient approximation (GGA) for exchange en¬ 
ergies. Electronic degrees of freedom were optimized with 
a conjugate gradient algorithm. Valence electron config¬ 
urations were: Mo_pv 4^^554/; Te_5^p^. Van der 

Waals interactions that bond the three-atom thick 2D X- 
Mo-X layers (X=S, Te) together were modeled with the 
non-local correlation functional of Klimes et al}^ To¬ 
tal energies were also calculated without van der Waals 
interactions, but, up to a basis of 140 structures, ground- 
state analyses always predicted false ground-states (sup¬ 
plementary material). Gonvergence with respect to k- 
point meshes was achieved by increasing the number of k- 
points until the total energy converged. A 500 eV cutoff- 
energy was used in the ’’high precision” option, which 
converges absolute energies to within a few meV/mol (a 
few tenths of a kJ/mol of exchangeable S- and Te-anions). 
Precision is at least an order of magnitude better. Resid¬ 
ual forces of order 0.02 eV or less were typical. Often, 
convergence with respect to hexagonal c-axis length was 
not automatic, and it was necessary to chose an initial 
c-axis value that is close to the converged value. Galcu- 
lated interlayer spacings in M 0 S 2 and MoTe 2 are 2.992 
A and 3.513 A, respectively, corresponding experimental 
values are: 2.977 A and 3.382 A 

Formation energies (AEf) for 233 Moi{SmP^n )2 su¬ 
percells are plotted in Fig. 1, in which values for AEf are 
normalized per mol of exchangeable anions, S and Te: 

AEf = {Estr - tuXmoS's “ tiFmoTcs)/( 2(7n + n)) (1) 
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FIG. 1. Comparison of formation energies,AF^/, for the 235 
DFT calculations (solid circles, green online) to Cluster Ex¬ 
pansion (CE) formation energies: AEnt (large open squares, 
red online) is the CE-fit to the DET set; AEgs (smaller open 
squares, blue online) are the CE-based ground-state analysis; 
AEi (solid black diamond dX X — 0.46, AEi ^ 0.03 eV) is 
the /-phase formation energy. All AE/>0 implies that there 
are no ordered ground-states, and suggests that the phase 
diagram will have a miscibility gap. 


Here: Estr is the total energy of the Mo/(SmTen )2 su¬ 
percell; EmoS 2 i^ energy/mol of M 0 S 2 ; ^MoTe 2 is the 
energy/mol of MoTe 2 . 

All supercell energies are positive which suggests a mis¬ 
cibility gap system, unless one or more entropy stabilized 
phases are stable. 


B. The Cluster Expansion Hamiltonian 

A cluster expansion Hamiltonian (CEH)^^, for the (1- 
X)-MoS 2 -(X)-MoTe 2 quasibinary system was fit to the 
set of 235 formation energies, AEvasp^ solid dots (green 
online) in Fig. 1 with a cross validation score of (CV 
)2=0.00723896). Fitting of the CEH was performed with 
the Alloy Theoretic Automated Toolkit (ATAT)^^’^^“^^ 
which automates most of the tasks associated with CEH 
construction. A complete description of the algorithms 
underlying the code can be found in^^. Large open 
squares in Fig. 1 (red online) indicate values of the 
235 AEpit that were calculated with the CEH. Smaller 
open squares {AEqs^ blue online) indicate the results of 
a ground-state analysis in which the CE was used to cal¬ 
culate formation energies for all ordered configurations 
with 16 or fewer anion sites, 151,023 structures. 


III. RESULTS AND DISCUSSION 

A first principles phase diagram (FPPD) calcula¬ 
tion was performed with grand-canonical, and canoni¬ 
cal, Monte Carlo (MC) simulations using the emc2 and 
phb codes which are part of the ATAT package^^”^^. 


Most phase boundaries were calculated with the phb 
program which uses equilibration tests to set the num¬ 
bers of equilibration- and MC-passes^^. To draw high- 
T extensions of the two-phase fields, and to locate the 
I ^ r transition, a 48x48x12 unit cell simulation box 
box was used, with 2000 equilibration passes and 2000 
MC-passes (see suplimentary material for comparisons 
of various equilibration- and MC-pass settings in calcu¬ 
lations of the I ^ r phase transition). The predicted 
phase diagram is shown in Fig. 2; where {M 0 S 2 ) denotes 
an MoS' 2 -rich solution phase, and similarly for {MoTe 2 ). 



FIG. 2. Calculated phase diagram. The isothermal line at 
82.5 K indicates that the /-phase, at A 0.46, is entropy 
stabilized (not a ground-state); Here: {M 0 S 2 ) and {MoTe 2 ) 
indicate MoAS' 2 -rich and MoTe 2 -rich solid solutions, respec¬ 
tively. Dotted lines in the region of the I' ^ disordered tran¬ 
sition indicate inferred extensions of calculated phase bound¬ 
aries. 


Kang et al.^ performed first principles phase diagram 
calculations for four dichalcogenide monolayer systems: 
MoSe 2 {i-x)Te 2 x, WSe 2 {i-x)Te 2 x, MoS 2 {i-x)Te 2 x and 
WS 2 {i-x)Te 2 x‘ van der Waals interactions, and fitting 
Their CEs wer fit to about 40 structures per system, 
and van der Waals interactions (with substrate) omitted. 
All systems were predicted to have miscibility gaps, and 
surprisingly, all consolute points are on the Te-rich sides. 
One expects the consolute point to be on the S-rich side, 
because it typically requires less energy to substitute a 
smaller S-ion into a larger Te-ion site, than vice versa. 
Figure 2 also has reduced solubility on the Te-rich side, 
but this is related to immiscibility between the I- (/')- 
phase, and the Te-rich phase. 

It is not clear that the Kang et al. M 0 S 2 — 
MoTe 2 phase diagram would still be a simple misci¬ 
bility gap had they included hundreds of structures in 
their CE-fit rather than about 40. Hence the monolayer 
vs. bulk comparison is uncertain. 

Surprisingly, the phase diagram predicted here has 
multiple two-phase fields, separated by two ordered in¬ 
commensurate phases, neither of which is a ground state. 
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To seven digits, the calculated bulk composition of the 
I-phase, just above its 82.5 K minimum temperature of 
stability, is X=0.4642857=13/28; i.e. Moi^Si^Tei^. Sta¬ 
bility of the I-phases is a robust result, in Monte-Carlo 
simulations: (1) CEH fits to 128, 153, 162, 182, 225 and 
235 formation energies all predict I-type ordering; (2) 
I-phase forms spontaneously on cooling of the /'-phase 
at T ^ 500//, and on heating of a low-T equilibrium 
{M 0 S 2 ) + {MoTe 2 ) assemblage to T ^ 82//; (3) I'- 
phase forms spontaneously on heating of the I-phase, or 
cooling of a disordered solid solution with X ^ 0.46. This 
calculation considers only Scon , and ignores excess vibra¬ 
tional entropy, Syib, which could conceivably destabilize 
the I-phases. In light of (1) above, however, this seems 
highly unlikely. Also, there is no fully satisfactory way to 
model Syib for an aperiodic phase, and a reasonable ap¬ 
proximate structure (with I-phase like ordering) would 
require at least a low symmetry 84-atom cell; which is 
beyond the scope of this study. 



FIG. 3. Minimum energies for various Monte-Carlo super¬ 
cells; is the length of the supercell, in c-axis units of the 
PGsImmc disordered-phase cell constants. The flat mini¬ 
mum at AEf = AEi ^ 0.03 eV/anion is interpreted as the 
/-phase formation energy. 


Figure 3 shows how the CE-calculated /-phase forma¬ 
tion energy AEf{I — phase) = AE/, varies as a func¬ 
tion of MC-supercell size and shape. The flat minimum 
at AEj ^ 0.03 eV/anion is the calculated /-phase for¬ 
mation energy which is plotted as the solid black dia¬ 
mond in Eig. 1. Supercell dimensions were chosen to 
accommodate Moi45'i5Tei3 stoichiometry and I-phase 
ordering. Note that many of the AEf plotted in Eigs. 
1, are lower in energy than AE/, but that they are for 
periodic structures with 16 or fewer anion sites in which 
Scon ^0 as T ^ 0//, and clearly (Eigs. 4) the /- 
phase is incommensurate, i.e. aperiodic with Scon > 
0. Eigures 4 exhibit the S:Te (yellowibrown online, re¬ 
spectively) ordering at: (a) 200 K; and (b) 575 K, i.e. 
below and above the / ^ /' phase transition. Mo-atoms 
are omitted for clarity, and labels (001)/:), (100)//, and 
(010)// refer to corresponding crystallographic planes in 
the high-T POs/mmc disordered-phase. 



FIG. 4. Monte-Carlo-snapshots of S:Te-ordering in the 
I- and /'-phases at: (a) X ^ 0.46 and T=200 K; 
and (b) X ^ 0.475 and T=575 K, respectively (online 
S=yellow, Te=brown, Mo omitted for clarity). Labels (001)//, 
(100)// and (010)//, refer to corresponding crystallographic 
planes in the high-T POs/mmc disordered phase. 


In the (001)/:)- and (100)//-planes, ...P^Ten... chains 
in the <010>// direction, most often have ra = 
5 or 6 and n = 4 or 5. Also, in (100)//, the ...SmTcn--- 
chains exhibit irregular alignments relative to one an¬ 
other. Note however, that P^Ten-chains in (001)//- 
planes order along the <001 >// direction, such that Po¬ 
unds alternate with Ton-units in adjacent 3-atom thick 
2D-layers; i.e. ...P^Te^P^/Te^/...chains are stacked on 
top of ...Pe^PmEon'Pm'•••chains with inescapable misfits, 
owing to the different and variable values of m and n. 
Thus /-phase ordering is inevitably imperfect, aperiodic, 
and incommensurate, which suggests that the I ^ I' 
phase transition is first-order. 

The difference between /- and /'-phases appears to be 
a distinction between 3D-ordering in the low-T /-phase 
and 2D-ordering in the high-T /'-phase. Clearly, ordering 
is stronger in the basal (001)//-plane than in the (100)//- 
or (010)//-planes, and striped order within (001)// per¬ 
sists above the I ^ I' transition. 

Figure 5a is a Monte-Carlo T-scan (heating) of the 
total energy Etot{T) which confirms first-order char¬ 
acter for the I ^ r transition: a critical (continu¬ 
ous) transition^^ would not exhibit the sharp change 
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FIG. 5. a) Monte-Carlo T-scans (heating) of: total energy, Etot{T), and Helmholtz energy, F{T), (inset) as functions of 
temperature; b) idealized schematic comparing equilibrium- and metastable transition paths (solid black line and red-dotted 
blue line, online, respectively). Etot{T) indicates a phase transition at T 500X; the absence of a clear change of slope in 
E(T) indicates that the transition is only weakly first-order; hence the dotted I ^ I' transition line in Fig. 2. 


at T ^ 503iF; also a transition from the lower-T, 
higher-phase to the higher-T, lower-E ’tot phase 
requires that the the lower-T phase be superheated, 
i.e. metastahle before it transforms. Figure 5b is 
an idealized schematic that compares equilibrium- and 
metastable-transition paths (solid black line, and blue 
line with red dots, online respectively). This transition 
is more subtle in cooling simulations, but still evident in 
snapshots. In the Fig. 5 inset, no change of slope at the 
transition is evident in the Helmholtz energy, F(T), which 
suggests that the transition is first-order^^. Also, 

when the MC temperature-increment was decreased from 
1.0 K/MC-step to 0.1 K/MC-step (not shown), the pre¬ 
dicted transition temperature decreased by about 10 K, 
which indicates more superheating at 1.0 K/MC-step 
than at 0.1 K/MC-step, hence first-order character. This 
transition is shown as a dotted line in Fig. 2 because the 
two-phase fields that a first-order transition implies are 
too narrow to resolve in the MC-simulations. 


IV. CONCLUSIONS 

To summarize, a first principles phase diagram calcula¬ 
tion for the 3D bulk system (1 — X) • M 0 S 2 ~ (-^) • MoTe 2 , 
that includes van der Waals interactions, predicts the for¬ 
mation of two entropy stabilized incommensurate^ i.e. 
aperiodic phases: the I- and /'-phases, at 0.46. 
Above the minimum temperature for stability of the I- 
phase, T ^ 82/F, the calculation predicts broad two- 
phase fields between the /-or /'-phase and disordered 
S- or Te-rich solution phases, {M 0 S 2 ) and (MoTe 2 ), re¬ 
spectively. Both the I ^ r and /' ^ disordered tran¬ 
sitions are predicted to be first-order. Dramatic changes 


in phase relations can be induced by arbitrarily small 
differences in energy, hence van der Waals interactions 
should not be ignored in layered 2D-systems such as 
M 0 S 2 — MoTe 2 . 
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SUPPLEMENTARY MATERIAL 


Mo(S,Te)- 



FIG. 1. Ground state analyses for cluster expansions based on fits to 102, 110, 120 and 140 formation energies, that were 
calculated without including van der Waals interactions. Note the presence of false ground-states in each fit; i.e. blue squares 
below the zero-energy line. 
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MoS,-MoTe, 



FIG. 2. Calculated I ^ I' transition points (solid up-pointing triangles, red online). Scatter in these results reflect MC- 
simulation fluctuations. The phase boundaries, {M 0 S 2 ) + ^ ^ + (MoTe 2 ), were calculated with 

the phb program in the ATAT package. 




X=0.464 
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T(K) 


FIG. 3. Calculated I ^ I' transition Etot(T) curves as functions of the numbers of MC- and equilibration-passes (n and eq, 
respectively) with a 48x48x12 MC-simulation box (solid circles), and a 56x56x12 MC-box (solid diamonds). 








